%solving strategy using residual based incremental 
%update static scheme with combination with 
%residual based newton raphson strategy
%this scheme supports for crack problems to enable controlling
%crack of element. One element is allowed to crack after each loop.
%in this scheme, lambda is changed during simulation to ensure maximum
%sigma1 is ft
function [elements,output] = solving_strategy_static_4...
    (elements,model,lambda0,dlambda,numTimeStep,maxits,...
    tol,postProcessInfo)
% tol = 1e-6;
% maxits = 50;
dim = 2;
output = [];
format long;
processInfo.integrationType = 'implicit';
processInfo.integrationStage = 'implicit';
disp('perform adaptive implicit Newton-Raphson integration scheme');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time steps iteration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
model.preAnalysis();
mesh = model.getMesh();
nDofs = mesh.countDofs();
r = zeros(nDofs,1);
r = model.assembleR(r);
u = zeros(nDofs,1);
lambda = lambda0;
lambda_old = 0;
sigma_max_old = zeros(model.countElements(),1);
for i=1:numTimeStep
    disp(['time step ' num2str(i) ' starts with lambda = ' num2str(lambda)]);
    lr = lambda*r;
    %update model's states
    model.setX(lambda,u);
    [elements,processInfo] = InitializeSolutionStep...
        (elements,model,processInfo);
    
    %assemble ri
    ri = assembleRi(elements,model);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Newton-Raphson loop
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    cnt = 0;
    %save the elements state
    elements_old = elements;
    %save the displacement
    u_old = u;
    while 1
        cnt = cnt+1;
        disp(['substep ' num2str(cnt) '--------------']);
        %assemble kt
        kt = assembleKt(elements,model);
%         rcond(kt)
%         find(eig(kt) < 0)
        
        du = kt\(lr-ri);
        u = u + du;
        model.setX(lambda,u);
        
        %update model's states
        [elements,processInfo] = CalculateMaterialResponse...
            (elements,model,processInfo);
        
        %assemble ri
        ri = assembleRi(elements,model);

        %compute convergent criteria
        %res = norm(du)/norm(u);
        res = norm(r-ri);
        
        disp(['substep ' num2str(cnt) ' complete']);
        if res < tol
            disp(['Newton Raphson iteration converged in ' ...
                num2str(cnt) ' loops']);
            %adaptive part
            [imax,max,list,maxarray] = findMaximumPrincipalStress(elements);
            if max > elements(imax).ft*(1+2e-6)
                %modify lambda
                %lambda = lambda*(elements(imax).ft*(1+1e-6)/max);
                lambda = lambda_old + ...
                    ((elements(imax).ft*(1+1e-6)-sigma_max_old(imax))/...
                    (max-sigma_max_old(imax)))*...
                    (lambda - lambda_old);
                disp(['adjusting lambda to ensure one element exceed ft, '...
                    'lambda_new = ' num2str(lambda)]);
                %restore the elements state
                elements = elements_old;
                %restore u
                u = u_old;
                %reset cnt
                cnt = 0;
            else
                sigma_max_old = maxarray;
                lambda_old = lambda;
                disp(['sigma_max_old = ' num2str(sigma_max_old(list))...
                    ', lambda_old = ' num2str(lambda_old)]);
                break;
            end
        end
        if cnt>maxits
            disp(['Newton Raphson iteration did not converge in ' ...
                num2str(maxits) ' loops']);
            break;
        end
    end
    
    [elements,processInfo] = FinalizeSolutionStep...
        (elements,model,processInfo);
    
    output = processOutput(i,dim,elements,model,...
        output,postProcessInfo,processInfo);
    
    disp(['time step ' num2str(i) ...
        ' finished with lambda = ' num2str(lambda)]);
    lambda = lambda + dlambda;
    disp('--------------------------------------------------------');
end

end
